Relaxation dynamics of the Kondo lattice model 
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We study the relaxation properties of the Kondo lattice model using the nonequilibrium dynamical mean 
field formalism in combination with the non-crossing approximation. The system is driven out of equilibrium 
either by a magnetic field pulse which perturbs the local singlets, or by a sudden quench of the Kondo coupling. 
For relaxation processes close to thermal equilibrium (after a weak perturbation), the relaxation time increases 
substantially as one crosses from the local moment regime into the heavy Fermi liquid. A strong perturbation, 
which injects a large amount of energy, can rapidly transform the heavy Fermi liquid into a local moment state. 
Upon cooling, the heavy Fermi liquid reappears in a two-stage relaxation, where the first step opens the Kondo 
gap and the second step corresponds to a slow approach of the equilibrium state via a nonthermal pathway. 

PACS numbers: 71.10.Fd, 71.27.+a 



I. INTRODUCTION 

Heavy Fermion compounds contain strongly interacting 
/-electrons which hybridize with extended s-, p- and d- 
electrons. The /-electrons are in a well-defined charge state, 
and at high temperature they act as local magnetic moments 
which scatter the conduction electrons. At low temperature, 
the hybridization between /- and conduction electrons can 
lead to the emergence of nontrivial electronic phases, such as 
the Kondo insulating state at half filling, or a strongly renor- 
malized Fermi liquid in the doped caseJ^ 

A simple model which captures essential aspects of the 
physics of these materials is the Kondo lattice model. 3 In this 
model, charge fluctuations of the /-electrons are completely 
suppressed. The localized degrees of freedom are described 
by spins S = 1/2, which are coupled to the spin of the con- 
duction electrons at the same site via an exchange interac- 
tion J (the Kondo coupling). For antiferromagnetic coupling 
(J > 0), a large value of J favors the formation of singlets 
on each site. At half-filling, this leads to the opening of a 
(pseudo-)gap in the spectral function of the conduction elec- 
trons. In the ferromagnetic Kondo lattice model (J < 0), a 
metallic phase is realized at small coupling, with a phase tran- 
sition (at half filling) to an insulating state at a critical value 
of J A 

The most interesting behavior is found away from half- 
filling, where an antiferromagnetic interaction with the local- 
ized moments leads to a renormalized bandstructure at low 
temperature, which resembles a flat band hybridized with the 
wide band of the conduction electrons. This Fermi liquid is 
characterized by strong mass renormalizations and a "large" 
Fermi surface, i.e., both conduction electrons and localized 
moments participate in the formation of the Fermi liquid state, 
and the Luttinger volume thus contains the total number of 
c- and /-electrons.- As the temperature is raised or the cou- 
pling J is reduced, a crossover occurs to a metallic state with a 
blurred, "small" Fermi surface, whose Luttinger volume con- 
tains the conduction electrons only. This physics has been 
beautifully demonstrated in a recent series of papers by Otsuki 
and collaborators, based on the dynamical mean field approx- 



imation (DMFT)££ 

In the present paper, we use the nonequilibrium extension 
of DMFT to investigate the real-time dynamics of the Kondo 
lattice model under strong nonequilibrium conditions. In par- 
ticular, we are interested in the timescales on which the sys- 
tem can undergo a transition between states with a large and 
small Fermi surface. In practice, it is easy to perturb the sys- 
tem so strongly that the heavy Fermi liquid is destroyed after 
re-thermalization at higher energy. We will implement such 
a perturbation by a sudden change of the interaction J, or a 
short magnetic field pulse, and investigate the crossover into 
the state with small Fermi surface in real time. The transi- 
tion back to the heavy Fermi liquid is then achieved upon 
cooling, by coupling the system to a dissipative environment 
which absorbs the energy injected into the system by the per- 
turbation. Although the maximal cooling rate is limited by the 
coupling strength, the way in which the Fermi liquid state is 
approached and the timescale for this process can still reveal 
intrinsic properties of the Kondo lattice model. To understand 
those relaxation times, we also consider the relaxation of the 
system close to equilibrium, by perturbing the local singlets 
with a weak magnetic field pulse. Due to the small amount of 
energy injected by such a weak pulse, the time evolution takes 
place only within one given phase and allows us to extract 
the equilibrium relaxation rate in the various temperature and 
doping regimes, and to connect this quantity to the presence 
or absence of strongly renormalized quasi-particles. 



II. MODEL AND METHOD 

The spin-i Kondo lattice model describes conduction elec- 
trons c interacting with localized electrons /. If the /-orbitals 
are half-filled and fluctuations into empty or doubly occupied 
states are energetically very expensive, only the spin degree of 
freedom S = ^ip^crifif remains in a low-energy model, while 

the charge fluctuations are suppressed [ijj^ = (ft, ft}]. The 
Hamiltonian of the Kondo lattice model then becomes 



2T = -5> 



$^ra«,<7 + J^Si-Si. (1) 
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FIG. 1: Top panel: illustration of the Kondo lattice model, describing 
conduction band electrons c hopping with matrix element v between 
orbitals (circles) and interacting via an exchange coupling J with the 
spin of the localized /-electrons (arrows). Bottom panel: dynami- 
cal mean field approximation of the Kondo lattice model, consisting 
of one c-electron orbital and the associated /-electron spin. The c- 
electrons couple to a self-consistently determined bath of noninter- 
acting sites, with hybridization function A(t, t'). 

Here, the first term corresponds to the kinetic energy of the 
conduction electrons, the second term gives the chemical 
potential contribution of the conduction electrons (n i a = 
cJo-C^ct), and the last term describes the interaction of the 
spins Si = \i\y c iCnpc^ of the conduction electrons with the 

localized electrons via the Kondo coupling J [ipl = (c|, cj)]. 
In this study, we will restrict our attention to antiferromagnetic 
J and to paramagnetic solutions. 

To investigate the properties of this model, we use dynam- 
ical mean field theory (DMFT) J. This approximate method, 
which becomes exact in the limit of infinite coordination 
number,- maps the lattice model onto a self-consistent so- 
lution of a quantum impurity model. The mapping can be 
applied to nonequilibrium situations 10,11 - by reformulating the 
theory on the Keldysh time contour C. As illustrated in Fig. Q] 
the impurity consists of a c-electron orbital coupled to a spin 
S and a bath of non-interacting sites. The relevant properties 
of the bath are encoded in the hybridization function A(t, t'), 
which describes the probability for transitions of c-electrons 
from the impurity site into the bath and back. Hence, the im- 
purity action reads 

S = -i JdtH l0C (t)-i Jdtdt'J2cUt)A(t,t') CtT (n ( 2 ) 

where 

H ioc = ^(n t + ni) + JS ■ -^l<T%j) c (3) 

is the local part of the Hamiltonian (Q3, consisting of the con- 
duction electron orbital and the spin S. 



We will consider a lattice whose noninteracting density of 
states is semi-elliptical with bandwidth 4u. In this case, the 
DMFT self-consistency becomes 

A(t,t')=v 2 G c (t,t'), (4) 

where G c (t,t') denotes the c-electron Green function (see 
Appendix). We set v = 1 as the unit of energy, and mea- 
sure time in units of u — 1 . The Green function G c (t,t') = 
—j(Tcc(t)c^ (<:')) must be computed numerically from the 
contour-ordered average (...) = Tr[Tce' s ...]/Tr[Xce' s ], using 
the action defined in Eq. (|2]i (Tc is the time-ordering operator 
on the Keldysh contour C). 

For equilibrium calculations, several techniques are avail- 
able to solve the impurity problem, among them the numeri- 
cal renormalization groupie and continuous-time Monte Carlo 
methods. — The latter come in two variants, the weak-coupling 
CT-J solver— and the hybridization expansion approach (CT- 
HYB^ In a CT-J simulation, the partition function of the im- 
purity model is expanded in powers of J, so that the Monte 
Carlo configurations consist of arbitrary sequences of spin- 
flip processes, with weight proportional to a determinant of a 
matrix of bath Green functions. In a CT-HYB calculation, the 
local problem iii oc is solved exactly, while the expansion of 
the partition function is done in powers of the hybridization 
function A. Here, the weight is proportional to the determi- 
nant of a matrix of hybridization functions. 

Neither CT-J nor CT-HYB suffers from a sign problem 
in equilibrium calculations. Monte Carlo simulations on the 
real-time Keldysh contour, however, lead to a dynamical sign 
problemr^i^ which restricts the simulations to rather short 
times. Since we are interested in both the transient dynam- 
ics of the Kondo lattice model and the long-time relaxation 
towards an equilibrium state, we use the non-crossing approx- 
imation (NCA) as an impurity solver. Similar to CT-HYB, 
NCA is based on an expansion of the impurity partition func- 
tion in powers of the hybridization functions. But rather than 
combining various (crossing and non-crossing) diagrams into 
a determinant, only the non-crossing diagrams are retained 
and summed up analytically via a Dyson equation. The imple- 
mentation of the NCA impurity solver on the Keldysh contour 
has been explained in Ref.[l7|, and we use the same procedure 
with the local Hamiltonian H\ oc represented as a 8 x 8 block 
matrix in the basis \S; n^-; nj_), with S =f, I and n a = 0, 1. 
Other formulations of (extendend) NCA for the Kondo lat- 
tice model have been previously usedJ£ An advantage of the 
present formulation with an 8-dimensional local problem is 
that it captures the formation of local singlets at the 0th order 
of the approximation, 19 The NCA solver does not suffer from 
a sign problem, so that the computational effort grows poly- 
nomially (like the third power) with the maximum simulation 
time. We converge the DMFT equations on the real-time con- 
tour time-step by time-step, using the procedure detailed in 
RefM 

In order to characterize the various equilibrium and 
nonequilibrium phases, we compute both static observables, 
such as the momentum distribution, and the frequency depen- 
dent spectral function. For a general time-evolving state, the 
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FIG. 2: Spectral functions of the n c — 1 Kondo insulator at ft = 50 
and indicated values of J. 




FIG. 3: Comparison between Green functions from CT-HYB and 
NCA for J = 1.5, n c — 1, and indicated temperatures. 



latter can be defined as 



1 f°° 

A(u, t) = Im / dt'e^ -^G^Ht', t), 

k Jt 



(5) 



where G Kt (t, t') = -i0(t - t')({c(t), c(t')}) is the retarded 
Green function. In practice, the time-integral is limited to 
some maximal time t' = t max , which could lead to artificial 
oscillations in the results. Below, this effect is well controlled 
with a suitably large cutoff t m . dx 50. 

In principle, there is always some arbitrariness in the def- 
inition of a nonequilibrium spectral function, and in general, 
the Fourier transform Eq. (0 need not even be positive. How- 
ever, in particular when the inverse width of the spectral fea- 
tures is small compared to the timescale on which A(u>,i) 
changes, the function ([5]l is closely related to a time-resolved 
photoemission and inverse photoemission spectrum.- 1 More- 
over, A(u>, t) constitutes a complete representation of the local 
Green function, and for an equilibrium state A(uj) = A(u>, t) 
becomes time-independent and reduces to the conventional 
definition. 



HI. RESULTS 
A. Equilibrium properties 

In order to set the stage for the study of the relaxation dy- 
namics and to test the quality of the NCA approximation, we 
first compute various results for the Kondo lattice model in 
equilibrium. Figure [2] shows the conduction electron spectral 
function @ at half-filling (p, = 0, n c = 1), for several values 
of J and inverse temperature f3 = 50. For sufficiently low 
temperature, the singlet formation between the conduction- 
electron and localized spins leads to the opening of a gap 
(Kondo insulator). The separation between the peaks is given 
by E(n = 2) + E(n = 0) - 2E(n = 1) = 1.5J, where 
E(n = 0) = 0, E(n = 1) = -f J - \i and E(n = 2) = -2/x 
are the lowest energy states for the local problem (H\ oc ) with 
n c-electrons.— The side-peaks which are split off by J cor- 
respond to the insertion or removal of an electron with addi- 
tional singlet-triplet excitations. Due to the exponential decay 
of the hybridization function in the Kondo insulating phase, 
we expect the NCA approximation to be rather accurate in 
this regime As shown in Fig. [3] already for J = 1.5 the 
NCA solution indeed provides a good approximation of the 
exact Green function, although NCA slightly underestimates 
the size of the Kondo gap, in contrast to the Mott gap in the 
Hubbard models 

As one dopes the system, a narrow quasi-particle peak ap- 
pears near the Fermi level at low temperatures. Figure JJshows 
spectral functions for J = 1.5, n c = 1.1 and n c = 1.4, 
and different inverse temperatures (3. As the temperature is 
increased, the narrow feature disappears ((3 = 10), but the 
Fermi level remains at the upper edge of the partially filled-in 
Kondo insulator gap. Eventually, the upper gap edge moves 
away from the Fermi level (/3 = 5) and at even higher tem- 
peratures (/3 = 2), the gap starts to fill in. At larger doping, 
the pseudo-gap is less pronounced and the low-temperature 
quasi-particle peak is merged with the upper band. 

To get a better understanding of the crossover between the 
various phases, we plot in Fig.[5]the temperature-dependence 
of the real part of the conduction electron self-energy at fre- 
quency to = 0, ReS(0), and the occupation p s i ng iet of the 
impurity singlet state (p s ingiet = {Pi(j ~ s ■ S)), where 
Pi = + — 2n^ni is the projector on the one-particle 
sector of the local Hilbert space). The behavior of ReS(0) 
was used in Ref. to define the crossover scale T* below 
which a coherent Fermi liquid state is formed. The evolu- 
tion of ReE(0) in our case looks similar to what was found 
in Ref. |^ for a hypercubic lattice, and the temperature depen- 
dence is also qualitatively consistent with the numerically ex- 
act CT-HYB results for the semi-circular DOS. We have used 
a linear extrapolation procedure to estimate ReE(0) from the 
values Yj(iu! n ) at non-zero Matsubara frequencies, rather than 
the quadratic fit employed in Ref.la, because this seems more 
appropriate at the relatively high temperatures considered in 
this study. At the filling n c = 1.1, the heavy Fermi liquid ap- 
pears below T* w 0.1, so that we can associate the shift of the 
Fermi level in Fig.|4]from a position within the (pseudo-)gap 
into the upper band and the formation of a narrow resonance 
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FIG. 4: J — 1.5. Temperature-dependence of the spectral function 
for n c = 1.1 (top panel) and n c = 1.4 (bottom) panel. 



with the appearance of heavy quasi-particles and the forma- 
tion of the large Fermi surface. Comparison with the tem- 
perature dependence of p s i n giet shows that the enhanced sin- 
glet formation sets in already at a higher temperature, below a 
crossover scale Tk ~ 0.5. Hence there is a temperature range 
T* < T < Tk between the Fermi liquid and local moment 
regime, where the singlet formation leads to a pronounced 
pseudo-gap in the spectral function and the Fermi level lies 
inside this gap. Above Tk, Psingiet decreases and the pseu- 
dogap in the spectral function starts to fill in. The existence 
of the two temperature scales T* and Tk in the single-site 
DMFT solution of the Kondo lattice model has been previ- 
ously discussed in Ref. |23l The crossover temperatures are 
correctly reproduced by the NCA approximation and, at least 
in the range considered here, they show little dependence on 
doping. 

The most direct evidence for a large Fermi surface is ob- 
tained from the momentum distribution n(e), which is given 
by the expectation value (c\,Ck) for band energy e = tk (see 
Appendix). It is shown in Fig. [6] for different couplings J 
at f3 = 50 and n c = 1.1, and 1.4. At small values of J, a 
smeared-out step is visible around the location of the Fermi 
surface of a free conduction electron gas. The formation of 
the heavy quasi-particle band leads to a step in the distribu- 
tion function n(e) at a different location and with a smaller 
size, corresponding to a shift of the Fermi surface and a re- 




FIG. 5: J = 1.5. Temperature dependence of the extrapolated value 
ReE(O) (left axis), showing the crossover from the Kondo insulator 
regime to the Fermi liquid regime below T*. Also plotted is the 
temperature dependence of the occupation of the singlet state (right 
axis), which illustrates the crossover from the local moment regime 
at high T to the Kondo insulator regime below Tk ■ The top panel 
shows data for n c = 1.1 and the bottom panel for n c — 1.4. 



duction of the quasiparticle weight. 

We conclude from these results that our NCA approxima- 
tion, which exactly treats an 8-dimensional local problem, 
provides a qualitatively correct description of the Kondo in- 
sulator, the heavy Fermi liquid and local moment regimes in 
the doped model, and of the various crossover phenomena. 
In the following, we will use this approach to study the re- 
laxation properties of the Kondo lattice model in the different 
parameter regimes. 



B. Relaxation close to thermal equilibrium 

We next compute the relaxation times of the system after it 
is weakly perturbed in one of the various parameter regimes 
discussed above. This will later on be important to understand 
the behavior of the system after a strong perturbation, for 
times long enough that a new equilibrium state is approached. 
The relaxation times are still entirely determined by equilib- 
rium properties of the corresponding electronic phase, which 
is not left by the weak perturbation. Technically, the following 
investigation could thus be done by computing a suitable lin- 
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FIG. 6: Distribution function n(e) for the antiferromagnetic Kondo 
lattice model at f3 = 50. The top panel corresponds to a filling n c = 
1.1 and the bottom panel to n c = 1.4. For large J, one observes the 
formation of a small kink in n(e), corresponding to the large Fermi 
surface, where the heavy quasiparticle band crosses the Fermi energy. 
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FIG. 7: Relaxation of p S m g i c t after a /i-field pulse of strength h = 
0.1, as a function of temperature. The initial state is the equilibrium 
state for indicated temperatures, n c — 1.1, J = 1.5. After the 
pulse, the system with /3 = 50 relaxes to a thermal state with inverse 
temperature « 15, and the system with /3 = 5 to a thermal state with 
inverse temperature sa 4.8. 



ear response quantity within the Matsubara formalism, but it is 
more convenient to do an explicit real-time calculation, which 
avoids analytical continuations and also does not require the 
evaluation of vertex functions of the impurity model. 

Specifically, we consider the relaxation of the system after 
a short and weak magnetic field pulse, which acts only on the 
localized spin by means of an additional perturbation h x (t)S x 
in the Hamiltonian ([3). The pulse is of the form 

h x (t) = h sin 2 (ir t/t pa i se ) (6) 

for t < <p U i S e and h x (t) = for t > t pu ise- We use pulses 
of duration t pu \ se = 1.5 and small field strength h, such that 
the heating effect is comparatively weak and we can study the 
relaxation time within the three temperature regimes T < T* 
("Fermi liquid", FL), T* < T < T K ("Kondo insulator", 
KI) and T > T K ("local moment", LM). We will focus on 
the time evolution of p s i ng iet. Other observables, such as the 
double-occupancy on the c-site, or the distribution function 
n(e, t) seem to give consistent results for the relaxation time. 

Because the pulse is acting only on one of the two spins 
which form the singlet (the localized spin S), it perturbs the 
singlet and leads to an initial decrease inp s j ng i et . This decrease 
is followed by a complicated transient evolution up to about 



t « 10, after which the system eventually settles into an expo- 
nential relaxation towards a new thermal equilibrium state at 
somewhat higher temperature. To measure the relaxation time 
r, we fit this long-time behavior with an exponential function 
Psingiet(t) = Psin g iet(i = oo) + Aexp(-t/r). For several sets 
of parameters we have cross-checked that within the accuracy 
of our calculation the extrapolated final value p s ingiet(* = oo) 
obtained from the fit corresponds to the value computed inde- 
pendently by assuming thermalization at constant energy. 

FigureQshows the time-evolution of p s ingiet(0 — Psingiet(i = 
oo) after a pulse of strength h = 0.1 and duration i pu i se = 1.5, 
for J = 1.5 and n c = 1.1. The initial state is the equilibrium 
state at different temperatures. In the FL regime (T < 0.1), 
the transient is relatively smooth, and the exponential relax- 
ation towards the equilibrium state becomes remarkably slow. 
In the KI regime 1/7<T< 1/3, three things happen: (i) the 
transient exhibits a plateau with oscillations whose period is 
roughly proportional to 1/J, (ii) the exponential relaxation 
becomes substantially faster with increasing temperature, and 
(iii) the exponential relaxation sets in from a value of p s ingiet 
which is much closer to the thermal value than in the FL 
regime. Finally, in the LM regime, the relaxation is so fast that 
after the initial transient (which seems to take about t « 10, 
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FIG. 8: Relaxation time r after a /i-field pulse of strength 0.1, as a 
function of temperature. The initial state is the equilibrium state for 
the indicated temperatures and n c = 1.1, J = 1.5. For comparison, 
we also plot ReE(O) and p s i ng iet (with arbitrary rescaling). 




FIG. 9: Relaxation time r for J = 1.5 as a function of doping (left 
panel) and at half-filling as a function of J (right panel). The initial 
state is the equilibrium state at /3 = 50. The large error bars at small 
J are due to the appearance of oscillations on top of the exponential 
relaxation. 

independent of temperature), the system has already thermal- 
ized. Within the accuracy of our simulation, it then no longer 
makes sense to fit an exponential to the long-time evolution. 

The behavior in Fig. [7] is reminiscent of the relaxation 
dynamics found after an electric field pulse in the Hub- 
bard model, as one approaches the insulator-metal crossover 
regime from the Mott insulating side.^ 4 In the doped Kondo 
lattice model, however, the slow relaxation occurs in the heavy 
Fermi liquid regime, while the KI regime with a deep pseudo- 
gap is associated with faster relaxation. 

In fact, it is the disappearance of coherent quasi-particles, 
and not the filling-in of the pseuo-gap, which leads to the 
faster relaxation. To illustrate this, we plot in Fig. [8] the 
relaxation time as a function of temperature, and compare 
this curve to ReE(0) and p s ingiet- The relaxation time tracks 
ReS(0) (associated with the formation of the heavy Fermi liq- 
uid with a large Fermi surface) and not p s i ne iet (associated with 
the opening of the gap). Furthermore, the relaxation time de- 
creases if the Fermi liquid is weakened, for example by in- 
creasing the doping (left panel of Fig. |9). That increasing 



doping leads to a weakening of the heavy Fermi-liquid state is 
indicated by the evolution of ReS(0) shown in Fig. [5] which 
shows that the FL state is formed at lower temperature for 
n c = 1.4 than for n c = 1.1, and by the spectral functions 
in Fig. [4] which show that the Kondo gap gets filled in with 
increasing doping. In the exhaustion limit n c — > 2, the Fermi 
liquid coherence temperature T* should drop to zero^ 

At least in the case J = 1.5 considered here, the Kondo 
insulator relaxes much faster than the weakly doped FL (about 
the same relaxation time as in the weakly doped model in the 
KI regime above T*), For smaller couplings, the relaxation 
time of the insulator increases rapidly (right panel of Fig. [9), 
since in the small- J limit the slowest processes are ~ 1/J. 
In the limit of large J, the relaxation time grows because it 
becomes difficult to transform an excitation energy of order J 
into kinetic energy. 



C. Fast melting of the large Fermi surface 

In this section, we study the dynamical evolution of the 
heavy Fermi liquid with a large Fermi surface into a state with 
small Fermi surface. This process can easily be triggered in 
many ways, provided a sufficient amount of energy is injected 
into the system. We will suddenly reduce the Kondo coupling 
J, which both leads to a considerable heating of the system 
and weakens the local singlets. The two effects combined are 
expected to result in a crossover from the FL phase into the 
KI and LM regimes. Figure [10] shows contour plots of the 
momentum occupation n(e,t) after a quench from J = 1.5, 
/3 = 50, n c = 1.4 to J final = 1.25, 0.75 and 0.5. The small 
quench (Jfi na i = 1.25) leads to a minor shift and a thermal 
broadening of the large Fermi surface (indicated by the red 
contour lines). The slow drift of the contour lines indicates 
that the relaxation time is slow, as expected for relaxation pro- 
cesses within the FL regime. The quench to Jg na i = 0.75 
drives the system into the KI regime, where Fermi liquid co- 
herence is lost (with a corresponding shift in ReS(0)), but 
the Kondo gap is still present. The relaxation is faster, in 
accordance with the faster dynamics measured after a weak 
magnetic field pulse in this regime. Finally, the quench to 
jfinai = 0.5 brings the system into the LM regime, character- 
ized by a small Fermi surface. Melting of the large Fermi sur- 
face and steepening of the momentum distribution around the 
location of the small Fermi surface happen on the timescale 
of a few inverse hoppings, comparable to the relaxation time 
after weak perturbations within the LM phase. We observe 
no additional bottleneck connected to the destruction of the 
large Fermi surface when hopping and J are comparable in 
magnitude. 

The evolution of the system into the KI and LM phases is 
also evident from the behavior of p s ingiet and the shift of ReS. 
Despite the heating effect, p s i ng iet remains large after quenches 
to J = 1.25 and J = 1, and the crossover into the LM 
regime is evident in the form of a substantial drop of p s ingiet 
for J = 0.5 (lower left panel of Fig. [lOj. The shift of ReS, 
on the other hand is evident in the bottom right two panels 
(corresponding to the quenches to J = 0.75 and J = 0.5): 



7 












t=6 




t=i2 




„ ' > initial 




/JjrV, thermal - 






7/ 




■/ / 

y / 









FIG. 10: Top panels: Time evolution of the distribution function n(e, t) after a quench from J = 1.5, n c = 1.4, /3 = 50 to J = 1.25, 0.75 
and 0.5 (from left to right). The plots show contours of constant n, with red contours corresponding to values of e, which in the initial state 
are associated with the large Fermi surface discontinuity. Bottom left panel: Time evolution of psingiet. Middle panel: Time evolution of the 
spectral function after a quench to J = 0.75. After thermalization, the system is in the KI phase (J3 sa 12.5). Right panel: Time evolution of 
the spectral function after a quench to J = 0.5 (effective temperature (3 ~ 8.3, in the LM phase). 



A time-dependent shift of ReS(0) is equivalent to a shift of 
the chemical potential fi. In nonequilibrium DMFT calcula- 
tions, a sudden shift of fi by A/i results in a rigid shift of the 
spectral function © by — Afj, on the frequency axis. In the 
present case, the rapid decrease of ReE(0) leads to a shift of 
the quasi-particle peak by +A(ReS(0)) < 0. Thermalization 
of the spectral function is seen to occur approximately on the 
same time-scale as the momentum distribution function. 



D. Formation of the heavy Fermion state upon cooling 

Finally, let us consider the dynamical formation of the FL 
state out of the LM phase. In general, it is not possible to 
reach the FL phase from the LM phase by a sudden increase 
of J, since the strong heating caused by such a quench would 
destroy the Fermi liquid coherence. We will thus approach the 
problem in a different way, which is at the same time closer 
to possible experiments on heavy-Fermion materials. In such 
an experiment, one might rapidly destroy the FL phase by a 
strong excitation (see Sec. IIII CK and monitor its reappear- 
ance out of the excited LM phase while energy is dissipated 
from the electronic system to the environment. As long as 
the intrinsic thermalization times of the system are fast com- 
pared to this dissipation time, the system will evolve through 
a sequence of equilibrium states of decreasing temperature, at 
a rate that is set by the coupling to the environment. How- 
ever, when the thermalization time becomes long, the system 
will fall out of equilibrium, or experience a slowdown of the 
cooling dynamics. Since our investigation in Sec. IIII B I has 



demonstrated a large increase of the thermalization times in 
the low-tempeature FL phase, one might expect the cooling 
dynamics in the doped Kondo lattice to reveal such a non- 
trivial time evolution. In the following, we will excite the 
system with a strong magnetic field pulse ©, and look at the 
subsequent dynamics in the presence of an additional particle 
reservoir at fixed temperature. 

To model dissipation of energy to other degrees of freedom 
we couple a thermal particle reservoir with inverse tempera- 
ture f3 locally to each site of the lattice. Technically, this cor- 
responds to a change of the DMFT self-consistency condition 
from Eq. (@]i to 

A(t,t') = v 2 G c (t,t') + A p {t,t'), (7) 

where the hybridization function of the reservoir is given by 

A /3 (t,t') = -ij de 7 (e)[G e (t,t') - fpiefie^W, (8) 

with f 13(e) the Fermi function (see Appendix). A simi- 
lar set-up has been previously considered in studies of the 
nonequilibrium properties of the Falikov-Kimball^ and Hub- 
bard models^ For 7(e) we use a semi-elliptic function with 
bandwidth 16 and amplitude 7(0) = 7. 

In the following we will focus on values of 7 which are 
large enough to allow for a fast energy dissipation, but still 
so small that the equilibrium properties of the system remain 
qualitatively unchanged with respect to the isolated system. 
We will therefore first study the effect of the bath (with in- 
verse temperature /3 = 50) on the spectral function and mo- 
mentum distribution function in equilibrium. As shown in the 
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FIG. 1 1 : Spectral functions A{ui) (left panels) and distribution func- 
tions n(e) (right panels) for different couplings 7 to a heat bath with 
/3 = 50. Top panels: J = 1.5. Bottom panels: J = 2.5. 



left panels of Fig.QT] a stronger coupling 7 reduces the peak 
near u> = in the spectral function, which is associated with 
the heavy quasi-particle band. Spectral weight is added to the 
gap region. As a result, the step-like feature in the distribu- 
tion function n(e) (near e = — 1) gets smeared out and for 
J = 1.5, n c w 1.20 it is hardly evident anymore for 7 > 0.4 
(right panels). To obtain a more prominent large Fermi sur- 
face, but still keep a strong coupling 7 to the environment, we 
also consider data for J = 2.5 and n c « 1.26 (bottom pan- 
els). While the average density is affected by the coupling to 
the bath, this effect is comparatively small: n c increases only 
by about one percent for a coupling strength 7 = 0.8. 

The time evolution of the system after an intense magnetic 
field pulse with strength h x = 0.5 and duration t max = 1.5 is 
plotted in Fig.Q~2] The left panels show the weight of the aver- 
age occupation of the singlet state, and the right panels show 
the density. The larger the coupling 7, the faster these observ- 
ables relax back to approximately the thermal value: If one de- 
fines a relaxation timescale To. 002 for this initial fast relaxation 
as the time by which p s i ng iet reaches its thermalized value up 
to within ±0.002, one finds that these relaxation times To. 002 
scale approximately linearly with I/7 ( J = 1.5: To. 002 = 62, 
36, 26.5, 22 for 7 = 0.2, 0.4, 0.6, 0.8, respectively; J = 2.5: 
T0.002 = 48, 26, 19, 15.5 for 7 = 0.2, 0.4, 0.6, 0.8 ). 

The fast initial relaxation leads to an overshooting of the 
singlet occupation and density, and it is followed by a slower 
convergence to the true steady state. The fast dynamics can 
be identified with the formation of the Kondo gap (more pre- 
cisely, a pseudo-gap), while the slow relaxation is related to 
the appearance of the heavy quasi-particle band. To illustrate 
this fact, we plot in Fig. [13] the time-dependent spectral func- 
tion [Eq. (O] for several times after the pulsed For J = 1.5, 
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FIG. 12: Time-evolution of the probability of the singlet state (left 
panels) and the density (right panels) for indicated couplings 7 to 
the bath with /3 = 50. The perturbation is a magnetic field pulse of 
amplitude h x = 0.5 and duration t max = 1.5. A larger coupling 7 
leads to a faster relaxation back to the thermal value. The top panels 
show data for J = 1.5 and the bottom panels for J = 2.5. 



7 = 0.6 (top panels), first indications of a feature near uj = 
appear around t = 20 (which is the time needed for p s i ng iet to 
relax back to roughly the thermal value), with a well-formed 
peak and a fully established gap around t = 30. However, 
the comparison with the equilibrium spectral function (black 
curve) shows that the amplitude of this peak is enhanced com- 
pared to the equilibrium peak (Fig. [13] middle panels), and the 
recovery of the latter takes much longer (the thermal curve is 
not reached for t w 40, which is the largest time for which we 
can compute well-resolved spectra). The observed long times 
needed to restore the heavy quasi-particle band are consistent 
with the relaxation times t « 50 measured for weak pertur- 
bations of the FL (Fig. [HJ. 28 Together with the corresponding 
feature at the lower gap edge, the formation of this coherent 
quasiparticle peak constitutes the bottleneck in the relaxation 
process. 

The shape of the peak in A(uj, t) near u> = indicates that 
the system does not approach the heavy Fermi liquid through 
a sequence of equilibrium states with decreasing temperature. 
Instead, it first establishes some kind of fairly stable non- 
thermal "precursor" state, which slowly evolves into the heavy 
Fermi liquid. In fact, comparison of the time-dependent spec- 
tra to equilibrium spectra at higher temperatures (Fig. Qj] right 
panels) shows that for increased temperature the quasi-particle 
peak would be broadened (and thus its maximum is shifted to- 
wards higher frequency), while in the time-dependent spectra, 
the location of the peak stays almost constant for large times. 
The larger amplitude of the peak in the precursor state can also 
not be explained with the small, transient change in doping - 
this would rather suggest a less prominent quasi-particle peak, 
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since for large times the doping is slightly smaller than in the 
final state (Fig. \T% . To further illustrate the difference be- 
tween the non-thermal state with enhanced quasi-particle peak 
and the true FL equilibrium state, we plot in Fig. [I4]the time- 
evolution of the momentum distribution function n(ek, t). The 
application of the pulse destroys the step-like feature marking 
the large Fermi surface. The fast relaxation of the singlet oc- 
cupation and density back to the thermal values is associated 
with a partial recovery of the step in n(ek,t), but it is clear 
from Fig.[Hthat after t = 20 (7 = 0.6) or t = 50 (7 = 0.2) 
the distribution function has not yet thermalized. The thermal 
distribution function with its well-defined step feature is only 
recovered once the heavy quasi-particle band has been fully 
reconstructed. 



IV. SUMMARY AND CONCLUSION 

We studied the relaxation dynamics of the Kondo lattice 
model (restricted to paramagnetic phases) using the nonequi- 
librium dynamical mean field formalism and an NCA im- 
purity solver which exactly treats an eight-dimensional local 
problem consisting of a spin and its associated conduction 
electron orbital. This approach is well-suited to describe the 
formation of local singlets, which are favored by antiferro- 
magnetic J. Comparison to data obtained with the numer- 
ically exact CT-HYB solver showed that the NCA approxi- 
mation yields qualitatively correct results in equilibrium over 
a wide doping and temperature regime. In particular, it cap- 
tures the crossover from a local moment regime into the heavy 
Fermi liquid regime as temperature is lowered in the doped 
system, with two crossover scales Tk (opening of a pseudo- 



gap in the conduction electron spectral function) and T* (shift 
in the real part of the conduction electron self-energy and for- 
mation of a large Fermi surface). 

We determined the relaxation times in these various phases 
and crossover regimes by applying a weak magnetic field 
pulse which perturbs the local singlets. While these numbers 
are related to equilibrium properties of the system, they may 
still be nontrivial to extract from a conventional imaginary- 
time equilibrium calculation. In the doped system, the relax- 
ation time grows with decreasing temperature, approximately 
proportional to the real part of the conduction electron self- 
energy. The slow relaxation in the heavy Fermi liquid phase 
is thus clearly associated with the existence of coherent heavy 
quasi-particles, and not, for example, with the presence of a 
pseudo-gap in the conduction electron spectral function. The 
relaxation time was also found to depend strongly on doping, 
with long relaxation times in the weakly doped heavy Fermi 
liquid regime. In contrast to the Hubbard model, the relax- 
ation times in the (Kondo) insulating state are substantially 
shorter than in the weakly doped regime, at least for J com- 
parable to the hopping. 

To study the relaxation dynamics of strongly excited sys- 
tems we considered quenches of the Kondo coupling and 
strong magnetic field pulses. Such perturbations lead to a con- 
siderable heating and it is thus easy to simulate the destruction 
of the low-temperature heavy-Fermi liquid state. By comput- 
ing the time-dependent momentum distribution function for 
quenches from intermediate to small J, we could demonstrate 
the destruction of the step feature associated with the large 
Fermi surface within a time of a few inverse hopping and the 
shift from a large to a small Fermi surface on a slower time 
scale (the relaxation time of the thermal state). 
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FIG. 14: J = 2.5. Time-evolution of the momentum distribution 
function for couplings 7 = 0.2 (top) and 7 = 0.6 (bottom) to the 
bath with f) = 50. 



In order to demonstrate the formation of the large Fermi 
surface upon cooling, we simulated the time-evolution after a 
strong magnetic field pulse in the presence of a thermal par- 
ticle reservoir, which removed the excess energy injected by 
the pulse. While a strong coupling to the heat-bath leads to a 
faster relaxation, it also smears out the heavy quasi-particle 
band. Nevertheless, for large enough Kondo coupling and 
large enough doping, one can have a well-defined large Fermi 
surface at low temperatures and fast relaxation. By comput- 
ing the time-dependent momentum distribution functions and 
spectral functions we could show that the relaxation back to 
the thermal state happens in two stages: after a fast initial re- 
laxation (whose time-scale depends on the strength of the cou- 
pling to the bath), a precursor state to the heavy Fermi liquid 
is formed, as evidenced by the appearance of a peak near the 
Fermi energy in the conduction electron spectral function, and 
a partially reconstructed step feature in the momentum distri- 
bution function. This fast relaxation is followed by a slower 
dynamics (presumably on time scales controlled by the long 
relaxation time in the heavy Fermi liquid), which leads to the 
thermalization of the narrow quasi-particle band also in the 
close vicinity of the Fermi energy. 

The existence and the nature of the precursor state should 
certainly be corroborated and further studied in future investi- 
gations. For example, it would be desirable to systematically 
look at those states at a slightly weaker coupling to the en- 



vironment in order to reduce the influence of the bath on the 
FL properties, in particular the quasi-particle weight and FL 
relaxation rate. However, smaller dissipation requires larger 
simulation times, which are not accessible within our cur- 
rent (single-processor) implementation of the NCA equations. 
Furthermore, it would be desirable, although computationally 
quite expensive, to check the influence of the NCA approxi- 
mation through a comparison to real-time results from higher 
order implementations of the self-consistent strong-coupling 
formalism. 

The calculations in this paper were limited to the param- 
agnetic phases of the antiferromagnetic Kondo lattice model. 
Many heavy electron materials however are close to a mag- 
netic instability. It would be interesting to extend our study 
to symmetry broken phases in order to enable an interplay be- 
tween electronic and magnetic excitations. 
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Appendix A: DMFT self-consistency with bath 

In this Appendix we briefly explain how a thermal 
fermionic bath can be incorporated into the DMFT equations 
for a semielliptic density of states. 

The local c-electron Green function of the impurity model, 
G c (t,t') = -iTr[T c e s c{t)J(t')]/Tr[T c e s } implicitly defines 
the self-energy £ via the impurity Dyson equation 

G c (t,t') = [id t + »-A{t,t')--E(t,t')}- 1 . (Al) 

Here and in the following, time arguments are on the Keldysh 
contour C, and Tq is the time-ordering operator on C. We 
use the notation for Keldysh equations detailed in Ref. I20I . 
Momentum-dependent lattice Green functions G k (t,t') = 
— i(TcCfe(t)c^(f')) are then obtained from the lattice Dyson 
equation 

G fc (M') = [id t + fi-e k -Z(t,t')]-\ (A2) 

where e k is the noninteracting dispersion. From these func- 
tions the time-dependent momentum distribution can be ob- 
tained, 

n(e k ,t) = -iG<(t,t). (A3) 

For a lattice whose noninteracting density of states is semi- 
elliptical with bandwidth Av, the local (momentum aver- 
aged) lattice Greenfunction can be shown to satisfy the self- 
consistent equation 

G(t,t')=Y,G k (t,t') (A4) 
fe 

= [idt + n-v 2 G(t,t')-Y,{t,t')]- x . (A5) 



Note that the second equality does not depend on the form 
of the self-energy, but only on the distribution of the band- 
energies £fc.— Comparison with the impurity Dyson equation 
dAU then yields the standard DMFT self-consistency condi- 
tion, Eq. ©. 

If an additional fermionic reservoir is coupled to the lattice 
at every site, one has to modify the lattice Dyson equation 
iA2\ by adding the hybridization function Ap(t,t')Xo the free 
dispersion, 

G k {t,t') = [idt + ii- e k - Ap(t,t') - n^t')]- 1 . (A6) 
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Hence the closed form equation for the momentum averaged 
Green function becomes 

G(t,t') = [id t + ii-v 2 G{t,t')~ A^M')- W)r\ 

(A7) 

and, by comparison with Eq. (IA1 b . we obtain the DMFT self- 
consistency with fermionic bath, 

A(t,t') =v 2 G(t,t') + Ap(t,t'). (A8) 
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